test visualizazion of GreenLand masks

link <- "/data/Dagobah/greengrass/grassland_ger/01_CTM_GL_mask_17-19/GL_mask_2017-2019.tif"

mask_GL <- terra::rast(link)

plot(mask_GL)

test vis for tsi output

Import TSI images (single location)

# Full mosaic
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt.ovr"

#  base case
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.vrt.ovr"

# Single tile base case
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif"

tsi_base_raw <- rast(link)

# Single tile case day_16_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

tsi_day_16_sigma_8_16_32_raw <- rast(link)

# Single tile case day_8_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_8_sigma_8_16_32/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

tsi_day_8_sigma_8_16_32_raw <- rast(link)

# Single tile case day_16_sigma_4_8_16
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_4_8_16/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

tsi_day_16_sigma_4_8_16_raw <- rast(link)

tsi_list <- list(tsi_base_raw, tsi_day_16_sigma_8_16_32_raw, tsi_day_16_sigma_4_8_16_raw, tsi_day_8_sigma_8_16_32_raw)

tsi_obs_length <- list()
lapply(seq_along(tsi_list),
       function(i) 
       tsi_obs_length[i] <<-  length(names(tsi_list[[i]]))
       )

test vis of TSI

# Set a "plan" for how the code should run.
plan(multisession, workers = 4)

# This does run in parallel!
future_map(seq_along(tsi_list), function(i)  
  terra::plot(tsi_list[[i]][[10:10]] )
  )

lapply(seq_along(tsi_list),
       function(i) plot(tsi_list[[i]][[10:18]])
       )

data wrangeling of TSI images

Compute and set min and max values for tsi cases

lapply(seq_along(tsi_list),
       function(i) 
      tsi_list[[i]] <<-  setMinMax(tsi_list[[i]])
)
## [[1]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 765  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4076026, 4106026, 3134920, 3164920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif 
## names       : 19840~LND05, 19840~LND05, 19840~LND05, 19841~LND05, 19850~LND05, 19850~LND05, ... 
## min values  :       -6318,        1235,       -4413,        2248,       -1832,       -2821, ... 
## max values  :        8778,        8231,        8899,        7515,        7072,        8358, ... 
## 
## [[2]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 866  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4076026, 4106026, 3134920, 3164920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif 
## names       : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ... 
## min values  :      NaN,      NaN,      NaN,      NaN,    -6318,    -6318, ... 
## max values  :      NaN,      NaN,      NaN,      NaN,     8778,     8778, ... 
## 
## [[3]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 866  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4076026, 4106026, 3134920, 3164920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif 
## names       : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ... 
## min values  :      NaN,      NaN,      NaN,      NaN,      NaN,      NaN, ... 
## max values  :      NaN,      NaN,      NaN,      NaN,      NaN,      NaN, ... 
## 
## [[4]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 866  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4076026, 4106026, 3134920, 3164920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif 
## names       : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ... 
## min values  :      NaN,      NaN,      NaN,      NaN,    -6318,    -6318, ... 
## max values  :      NaN,      NaN,      NaN,      NaN,     8778,     8778, ...

remove all NA layers

# create filter out all layers that dont contain any values
tsi_na_index <- list()
lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_na_index[i] <<- minmax(tsi_list[[i]]) |> 
         as.data.frame() |> 
         slice(1) |> 
         pivot_longer(everything()) |> 
         mutate(id = row_number()) |> 
         na.omit()
)

# convert into index vector
lapply(seq_along(tsi_list),
       function(i) 
        tsi_na_index[[i]] <<- as.vector(tsi_na_index[[i]])
         
)

# subset tsi object with only layers containing values
lapply(seq_along(tsi_list),
       function(i) 
         tsi_list[[i]] <<- tsi_list[[i]][[ tsi_na_index[[i]] ]] 
)

    # filter out reflectance values below and above 0-100 % range
lapply(seq_along(tsi_list),
       function(i) 
         tsi_list[[i]] <- tsi_list[[i]] |> 
  filter(across(1:length(names(tsi_list[[i]])), ~ . > 0),
         across(1:length(names(tsi_list[[i]])), ~ . < 10000))
)

Check NA layer occurance

lapply(seq_along(tsi_list),
       function(i) 
       tsi_obs_length[[i]] - length(names(tsi_list[[i]]))
)

lapply(seq_along(tsi_list),
       function(i) 
print(paste( round(( tsi_obs_length[i] / length(names(tsi_list[[i]]))), digits=4), "% of layers consisted completely of NA values (n=", ( tsi_obs_length[i] - length(names(tsi_list[[i]]))), "out of",  tsi_obs_length[i], ")")) 

)

test for R tile merging

link_a <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

link_b <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0052/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

a <- rast(link_a)
b <- rast(link_b)

ab <- merge(a,b)

plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
     range=c(0,10000)
    )

object.size((ab))

writeRaster(ab, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif",
            # overwrite=TRUE,
            datatype = "INT4S",
            filetype='GTiff')


link_ab <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif"

ab <- rast(link_ab)
plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
     range=c(0,10000)
    )

plot raw images

#rasterVis::levelplot(tsi[[c(24)]], par.settings = viridisTheme)
#rasterVis::levelplot(tsi[[c(1,4,24,25:35)]], par.settings = viridisTheme)


lapply(seq_along(tsi_list),
       function(i) 
      plot(tsi_list[[i]][[c(1:10)]], legend=F, col = viridis(n=100,option="D"),
     range=c(0,10000)
    )
)

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL

gif image creation

for (i in c(1:100
            #length(names(tsi))
                     )) {
  
   # Step 1: Call the pdf command to start the plot
  png(file = paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/plot_", i, ".png"), width = 265, height = 125, units='mm', res = 100) #,   # The directory you want to save the file in
   # width = 12, # The width of the plot in inches
    #height = 12) # The height of the plot in inches
plot(tsi[[i]], legend=T, col = viridis(n=100,option="D"),
     range=c(0,10000),
    plg=list( title = paste(as.Date(names(tsi[[i]]), format="%Y%m%d")), title.cex=1.25)
)


# Step 3: Run dev.off() to create the file
dev.off()
}
knitr::include_graphics("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/ts_16_8_16_32.gif")

spectra visualization

data manipulation and subsample creation

tsi_df <- list()
lapply(seq_along(tsi_list),
       function(i) 
         tsi_df[[i]] <<- terra::as.data.frame(tsi_list[[i]][[c(1:100)]], xy = TRUE, na.rm=F) |>
         # filter out rows that only contain na
         filter_at(vars(-x,-y), any_vars(! is.na(.))) |> 
         # random subset for data vis
         slice_sample(n=10)
)

# create long version of df with unique and formatted location and date cols
tsi_df_long <- list()
lapply(seq_along(tsi_list),
       function(i) 
         tsi_df_long[[i]] <<- tsi_df[[i]] |> 
         pivot_longer(cols = c(3:ncol(tsi_df[[i]])), names_to = "timestamp") |> 
         mutate(xy= paste0(x,"_",y)) |> 
         mutate(date = as.Date(timestamp, format="%Y%m%d"),
                year = lubridate::year(date))
       
)

plot faceted time series

lapply(seq_along(tsi_list),
       function(i) 
ggplot(tsi_df_long[[i]], aes(date, value, group=xy, color=xy)) +
  geom_line() +
  geom_point(size=0.5) +
  theme_classic() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "1 month", date_labels = "%m")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="none") +
  facet_wrap(~year, scales = "free_x")

)
## [[1]]
## Warning: Removed 42 row(s) containing missing values (geom_path).
## Warning: Removed 434 rows containing missing values (geom_point).

## 
## [[2]]
## Warning: Removed 21 row(s) containing missing values (geom_path).
## Warning: Removed 177 rows containing missing values (geom_point).

## 
## [[3]]
## Warning: Removed 18 row(s) containing missing values (geom_path).
## Warning: Removed 259 rows containing missing values (geom_point).

## 
## [[4]]
## Warning: Removed 42 row(s) containing missing values (geom_path).
## Warning: Removed 247 rows containing missing values (geom_point).

plot time series with formatted dates

ggplot(tsi_df_long, aes(date, value, group=xy, color=xy)) +
  geom_line() +
  theme_classic() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "1 month", date_labels = "%m-%y")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="none")

ggplot(tsi_df_long, aes(date, value)) +
  geom_smooth() +
  theme_classic() +
    scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
  theme(axis.text.x = element_text(angle = 45 ,  hjust=1)) 

ggplot(tsi_df_long, aes(date, value, group = timestamp)) +
  geom_boxplot() +
  theme_classic() +
  scale_x_date( date_labels = "%d-%m-%y")+
  theme(axis.text.x = element_text(angle = 45, hjust=1))

ggplot + tidyterra vis option (nice but runs slow)

ggplot() +
  geom_spatraster(data = tsi[[c(4:8)]]) +
  facet_wrap(~lyr, ncol = 2) +
  scale_fill_whitebox_c(
    palette = "viridi",
   # labels = scales::label_number(suffix = "ยบ")
  ) +
  labs(fill = "NDVI") +
  theme_classic()